import java.math.BigDecimal;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;


public class Part1 
{
	public Part1()
	{
		
	}
	
	//TODO: PART 1
	/**
	 * Creates and returns a hilbert matrix of dimension n
	 * @param n The dimension of the hilbert matrix
	 * @return The hilbert matrix
	 */
	public double[][] createHilbertMatrix(int n)
	{
		double[][] hilbert = new double[n][n];
		for (int i = 0; i < n; i++)
			for (int j = 0; j < n; j++)
				hilbert[i][j] = ((double)1) / (i+1 + j+1 -1);
		
		return hilbert;
	}
	
	
	public void LUDecomp(double[][] h, double[] b)
	{
		
//		int vIndex = 0;
//		double pivot = 0;
//		//first step, have to create G matrix
//
//		for (int i = 0; i < h[0].length; i++)
//		{
//			pivot = h[i][i];
//
//			for (int j = i+1; j < h.length; j++)
//			{
//				double test = h[j][i]- h[j][i] * h[j][i]/pivot;
//				
//				if (i==0)
//					break;
//				if (i == j-1)
//					h[j][i] = h[j][i]- h[j][i]/pivot;
//				else
//					h[j][i] = h[j][i]- h[j][i] * h[j][i]/pivot;
//			}
//		}	
//		vIndex++;

		double[][] L = createIdentity(h.length);
		double[][] U = createIdentity(h.length);
		double[][] Q = createIdentity(h.length); 
		double[][] G = createIdentity(h.length);

		
		for (int a = 0; a < h.length; a++) //keep track of pivot
		{
			for (int i = a+1; i < h.length; i++)//go through all rows
			{
				for (int j = a; j < h[0].length; j++)//go through all columns
				{
					double element = h[i][j] - h[a][j] * h[a][i]/h[a][a];
					double pivot = h[a][a];
					//if (element <= Math.pow(10,-15))
					//	element = 1;
					double multiplicationFactor = h[a][j] * element/pivot;
					h[i][j] = element; //- h[a][j] * element;
				}
				
				G = createIdentity(h.length);
				G[i][a] = h[a][i]/h[a][a]; //inverse is to take the oppposite, normally this value is negative
				
				L = multiply(L,G);
			}
			
		}
		
		System.out.println();
		printMatrix(L);
		System.out.println();

		
		printMatrix(h);
		
		System.out.println("Original Matrix:");
		printMatrix(multiply(L,h));
		System.out.println("END");
	}
	
	//TODO: PART 2
	public void houseHolder(double[][] h, double[][] b)
	{
		double[][] herm = createIdentity(h.length);
		double[][] Q = createIdentity(h.length);
		double[][] temp = new double[h.length][h.length];
		double[][] R = h;
		
		for (int i = 0; i < R[0].length-1; i++)
		{
			temp = new double[R.length-i][R.length-i];
			//create submatrix to anaylze
			for (int j = i; j < R[0].length; j++)
			{
				for (int a = i; a < R.length; a++)
					temp[a-i][j-i] = R[a][j]; 
			}
			
			
//			System.out.println();
//			System.out.println("u:");
//			printMatrix(calculateU(temp));
//			System.out.println();
			
			temp = subtract(createIdentity(temp.length), scalarMultiply((multiply(calculateU(temp), transpose(calculateU(temp)))),2));
			
//			System.out.println();
//			System.out.println("temp:");
//			printMatrix(temp);
//			System.out.println();
			
			herm = createIdentity(h.length);
			
			for (int j = i; j < h[0].length; j++)	//h[0].length = width, h.length = height
			{
				for (int a = i; a < h.length; a++)
				{
					herm[j][a] = temp[j-i][a-i];
				}
			}
			
//			System.out.println();
//			System.out.println("herm:");
//			printMatrix(herm);
//			System.out.println();
			
			R = multiply(herm,R);
			Q = multiply(Q,herm);
			
//			System.out.println();
//			System.out.println("R:");
//			printMatrix(R);
//			System.out.println();
//			System.out.println();
//			System.out.println("Q:");
//			printMatrix(Q);
//			System.out.println();
		}
		
//		System.out.println();
//		printMatrix(multiply(Q,R));
//		System.out.println();
		
		double[][] resultantMatrix = multiply(transpose(Q),b);
		double[][] x = new double[R.length][1];
		double curTotal = 0;
		
		for (int i = R.length - 1; i >= 0; i--) //i keeps track of rows
		{
			curTotal = 0;
			for (int j = R[0].length -1; j >= 0; j--) //j keeps track of columns
				if (i == j)
					x[i][0] = (resultantMatrix[i][0] - curTotal)/ R[i][j];
				else if ( j < i)
					continue;
				else
					curTotal += R[i][j] * x[j][0];
		}

		System.out.println("SOLUTION HERE!");
		printMatrix(x);
		System.out.println("SOLUTION HERE!");
		
		System.out.println(maxNorm(subtract(multiply(Q,R),h)));
		System.out.println(maxNorm(subtract(multiply(h,x),b)));
		
	}
	
	/**
	 * Calculates the U matrix of the h matrix
	 * @param h
	 * @return
	 */
	public double[][] calculateU(double[][] h)
	{
		double[][] u = new double[h.length][1];
		
		for (int i = 0; i < h.length; i++)
				u[i][0] = h[i][0];
		
		if (h[0][0] > 0)
			u[0][0] = h[0][0] + norm(u);
		else if (h[0][0] < 0)
			u[0][0] = h[0][0] - norm(u);
		else if (h[0][0] == 0)
			u[0][0] = h[0][0];
		
		u = scalarMultiply(u,1/norm(u));
		
		return u;
	}
	
	//TODO: PART 3
	public void givens(double[][] h, double[][] b)
	{
		double[][] R = h;
		double[][] Q = createIdentity(h.length);
		
		//for (int a = 0; a < h.length; a++) //keeps track of the pivot row number
			for (int i = 0; i < h[0].length; i++) //look at each row, column stays constant
				for (int  j = i+1; j < h.length; j++)
				{
					//System.out.println();
					//printMatrix(calculateGivensG(R,i,j,i));
					//System.out.println();

					Q = multiply(Q, transpose(calculateGivensG(R,i,j,i)));
					R = multiply(calculateGivensG(R,i,j,i), R);
				}
			
//			printMatrix(R);
//			System.out.println();
//			printMatrix(Q);
//			System.out.println();
//			printMatrix(multiply(Q,R));
			
			double[][] resultantMatrix = multiply(transpose(Q),b);
			double[][] x = new double[R.length][1];
			double curTotal = 0;
			
			for (int i = R.length - 1; i >= 0; i--) //i keeps track of rows
			{
				curTotal = 0;
				for (int j = R[0].length -1; j >= 0; j--) //j keeps track of columns
					if (i == j)
						x[i][0] = (resultantMatrix[i][0] - curTotal)/ R[i][j];
					else if ( j < i)
						continue;
					else
						curTotal += R[i][j] * x[j][0];
						
				
			}
			
			System.out.println("SOLUTION HERE!");
			printMatrix(x);
			System.out.println("SOLUTION HERE!");
			
			System.out.println(maxNorm(subtract(multiply(Q,R),h)));
			System.out.println(maxNorm(subtract(multiply(h,x),b)));
	}
	
	/**
	 * Calculates the G matrix for matrix m for two elements
	 * @param m The actual matrix
	 * @param index1 The index of the first value to make a G matrix of
	 * @param index2 The index of the second value to make a G matrix of
	 * @return The G matrix rotation for this
	 */
	public double[][] calculateGivensG(double[][] m, int m1_row, int m2_row, int m_column)
	{
		double[][] g = createIdentity(m.length);
		
		for (int i = 0; i < g.length; i++) //traverses through columns, row number is constant ROW
			for (int j = 0; j <g[0].length; j++) //column
			{
				//STUFF NEEDS TO GO HERE, THE STUFF FOLLOWING THE IF STATEMENTS IS WRONG AND I WILL CORRECT THIS TODAY
				//3/30/11- 12:55 PM
				if (i == m1_row && j == m_column)
					//cos
					g[i][j] = m[m1_row][m_column] / (Math.sqrt(Math.pow(m[m1_row][m_column],2) + Math.pow(m[m2_row][m_column],2)));
				else if (i == m2_row && j == m_column)
					//sin
					g[i][j] = -1 * m[m2_row][m_column] / (Math.sqrt(Math.pow(m[m1_row][m_column],2) + Math.pow(m[m2_row][m_column],2)));
				else if (i == m1_row && j == m_column + m2_row-m1_row)
					//-sin
					g[i][j] = -1 * -1 * m[m2_row][m_column] / (Math.sqrt(Math.pow(m[m1_row][m_column],2) + Math.pow(m[m2_row][m_column],2)));
				else if (i == m2_row && j == m_column+m2_row-m1_row)
					//cos
					g[i][j] = m[m1_row][m_column] / (Math.sqrt(Math.pow(m[m1_row][m_column],2) + Math.pow(m[m2_row][m_column],2)));
			}
		
		return g;
		
	}
	
	//TODO: GENERAL USE METHODS
	/**
	 * Prints the matrix to the screen
	 * @param m The matrix to be printed to the screen
	 */
	public void printMatrix(double[][] m)
	{
		for (int i = 0; i < m.length; i++)
		{
			for (int j = 0; j<m[0].length; j++)
			{
				System.out.print(m[i][j]);
				System.out.print(" ");
			}
			System.out.println();
		}

	}
	
	/**
	 * Multiplies two matrices together and returns the result. 
	 * NOTE: The matrix {{3,2},{1,4}} is the following:
	 * 3 2 
	 * 1 4
	 * @param m1 The first matrix
	 * @param m2 The second matrix
	 * @return The resulting matrix
	 */
	public static double[][] multiply(double[][] m1, double[][] m2)
	{
		double[][] result = new double[m1.length][m2[0].length];
		
		  for(int i = 0; i < m1.length; i++) 
			  for(int j = 0; j < m2[0].length; j++) 
		    	  for(int k = 0; k < m1[0].length; k++)
		    		  result[i][j] += m1[i][k]*m2[k][j];
		  
		return result;
	}
	
//	/**
//	 * Inclusive indices
//	 * @param m
//	 * @param leftIndex
//	 * @param rightIndex
//	 * @param topIndex
//	 * @param bottomIndex
//	 * @return
//	 */
//	public double[][] subMatrix(double[][] m, int leftIndex, int topIndex, int rightIndex, int bottomIndex) {
//		
//		double[][] subMatrix = new double[bottomIndex - topIndex + 1][rightIndex - leftIndex + 1];
//		
//		for (int y = topIndex; y <= bottomIndex; y++) {
//			for (int x = leftIndex; x <= rightIndex; x++) {
//				subMatrix[y - leftIndex][x - topIndex] = m[y][x];
//			}
//		}
//		
//		return subMatrix;
//	}
	
//	private class rowReduceReturn {
//		double[][] reducedMatrix;
//		
//		List<double[][]> gMatrices;
//		
//		public rowReduceReturn(double[][] reducedMatrix,
//				List<double[][]> gMatrices) {
//			super();
//			this.reducedMatrix = reducedMatrix;
//			this.gMatrices = gMatrices;
//			if (this.gMatrices == null)
//				this.gMatrices = new ArrayList<double[][]>();
//		}
//	}
	
//	public double[][] getIdentityMatrix(int width, int height) {
//		double[][] m = new double[width][height];
//		
//		for (int i = 0; i < m.length; i++) {
//			Arrays.fill(m[i], 0);
//		}
//		
//		int i = 0;
//		while(i < width && i < height) {
//			m[i][i] = 1;
//			i++;
//		}
//		
//		return m;
//	}
	
//	public rowReduceReturn rowReduce(double[][] m) {
//		return rowReduce(m, m[0].length, m.length);
//	}
	
//	private rowReduceReturn rowReduce(double[][] m, int originalWidth, int originalHeight) {
//		List<double[][]> gMatrixList = new ArrayList<double[][]>();
//		
//		if (m.length <= 1)
//			return new rowReduceReturn(m, null);
//		
//		//double[][] gMatrix = getIdentityMatrix(m[0].length, m.length);
//		
//		//double factor = m[1][0] / m[0][0];
//		for (int y = 1; y < m.length; y++) {
//			double factor = m[y][0] / m[0][0];
//			for (int x = 0; x < m[0].length; x++) {
//				m[y][x] = m[y][x] - (m[0][x] * factor);
//				double[][] gMatrix = getIdentityMatrix(m[0].length, m.length);
//				gMatrix[y][0] = -factor;
//				gMatrixList.add(gMatrix);
//			}
//		}
//		
//		double[][] finalMatrix = Arrays.copyOf(m, m.length);
//		
//		rowReduceReturn subMatrixResult = rowReduce(subMatrix(m, 1, 1, m[0].length - 1, m.length - 1), originalWidth, originalHeight);
//		
//		//for each gMatrix, expand size by one
//		for (int i = 0; i < subMatrixResult.gMatrices.size(); i++) {
//			double[][] current = subMatrixResult.gMatrices.get(i);
//			double[][] temp = new double[current.length][current[0].length];
//			for (double[] t : temp)
//				Arrays.fill(t, 0);
//			temp[0][0] = 1;
//			for (int y = 0; y < current.length; y++) {
//				for (int x = 0; x < current[0].length; x++) {
//					temp[y+1][x+1] = current[y][x];
//				}
//			}
//			subMatrixResult.gMatrices.set(i, temp);
//		}
//		
//		for (int y = 0; y < subMatrixResult.reducedMatrix.length; y++) {
//			for (int x = 0; x < subMatrixResult.reducedMatrix[0].length; x++) {
//				finalMatrix[y+1][x+1] = subMatrixResult.reducedMatrix[y][x];
//			}
//		}
//		
//		gMatrixList.addAll(subMatrixResult.gMatrices);
//		return new rowReduceReturn(finalMatrix, gMatrixList);
//	}
	
	public double[][] scalarMultiply(double[][] m1, double val)
	{
		double[][] result = new double[m1.length][m1[0].length];
		
		  for(int i = 0; i < m1.length; i++) 
			  for(int j = 0; j < m1[0].length; j++) 
		    		 result[i][j] += m1[i][j]*val;
		  
		return result;
	}
	
	/**
	 * Creates and identity matrix of dimension n
	 * @param n The dimensions of the identity matrix
	 * @return The identity matrix
	 */
	public double[][] createIdentity(int n)
	{
		double[][] identity = new double[n][n];
		
		for (int i = 0; i < n; i++)
			for (int j = 0; j< n; j++)
				if (j == i)
					identity[i][j] = 1;
				else
					identity[i][j] = 0;
		
		return identity;
	}
	
	public double[][] createOnesMatrix(int n)
	{
		double[][] ones = new double[n][1];
		
		for (int i = 0; i < n; i++)
				ones[i][0] = 1;
		
		return ones;

	}
	
	public double[][] add(double[][] m1, double[][] m2)
	{
		double[][] add = new double[m1.length][m1.length];
		
		for (int i = 0; i < m1.length; i++)
			for (int j = 0; j< m1.length; j++)
				add[i][j] = m1[i][j] + m2[i][j];
		
		return add;
	}
	
	public double[][] subtract(double[][] m1, double[][] m2)
	{
		double[][] add = new double[m1.length][m1.length];
		
		for (int i = 0; i < m1.length; i++)
			for (int j = 0; j< m1[0].length; j++)
				add[i][j] = m1[i][j] - m2[i][j];
		
		return add;
	}
	
	public double[][] transpose(double[][] m1)
	{
		double[][] transpose = new double[m1[0].length][m1.length];
		
		for (int i = 0; i < m1[0].length; i++)
			for (int j = 0; j< m1.length; j++)
				transpose[i][j] = m1[j][i];
		
		return transpose;
	}
	
	public double norm(double[][] m1)
	{		
		double tempTotal = 0;
		
		for (int i = 0; i < m1.length; i++)
		{
			tempTotal+=Math.pow(m1[i][0],2);
		}
		
		return Math.sqrt(tempTotal);
	}
	
	/**
	 * Returns the maximum norm for error calculation
	 * @param m1 The matrix to calculate the max norm for
	 * @return The max norm of the matrix
	 */
	public double maxNorm(double[][] m1)
	{
		double maxNorm = Integer.MIN_VALUE;
		
		for (int i = 0; i < m1.length; i++)
			for (int j = 0; j <m1[0].length; j++)
				if (m1[i][j] > maxNorm)
					maxNorm = m1[i][j];
		
		return maxNorm;
	}
	
	//TODO: TESTING
	
	public static void main(String[] args)
	{
		Part1 test = new Part1();
		double[][] m1 = {{3,2},{1,4}};
		double[][] m2 = {{3,2},{1,4}};
		
		double[][] house_Test = {{1,3,3}
						, {0,3,2}
						, {2,0,1}};
		
		double[][] givens_Test = {{1,2,0}
		, {1,1,1}
		, {2,1,0}};
		
		double[][] generalQR_testA = {{2,1,1}
		, {4,-6,0}
		, {-2,7,2}};
		
		double[][] generalQR_testb = {{1},{1},{1}};
		
		test.printMatrix(test.createHilbertMatrix(5));
		test.printMatrix(test.multiply(m1, m2));
		
		test.createOnesMatrix(5);
		
		//test.houseHolder(test.createHilbertMatrix(5), null);
		test.houseHolder(test.createHilbertMatrix(5), test.createOnesMatrix(5));
		
		
		/*System.out.println();
		test.printMatrix(test.rowReduce(m1));
		
		double[][] m3 = {
				{ 1, 2, 3, 4, 5},
				{ 3, 4, 2, 5, 6},
				{ 5.4, 5, 6, 3.33, 2}, 
				{ 2, 3, 4, 1, 5 }
		};
		
		System.out.println();
		test.printMatrix(test.rowReduce(m3));*/
		
		double[][] m4 = {
				{ 2, 1, 1},
				{ 4, -6, 0},
				{ -2, 7, 2}
		};
		
		System.out.println();
		System.out.println();
		//test.printMatrix(test.rowReduce(m4).reducedMatrix);
		
		test.givens(generalQR_testA, generalQR_testb);
		test.houseHolder(generalQR_testA, generalQR_testb);

		
		System.out.println();
		
		//test.givens(test.createHilbertMatrix(20),test.createOnesMatrix(5));
		
		test.LUDecomp(test.createHilbertMatrix(5), null);
	}
}
